NASA Contractor Report 1 89065 


e^i 


Developments in Variational Methods 
for High Performance Plate 
and Shell Elements 


(NASA-CR-1&9065) OFVFLQPMENTS IN 
VAR I AT I FINAL MFTHQOS FOR HIGH PERFORMANCF 
PLATE AND SHtLL ELEMENTS Final Report, Jan. 
1989 (Colorado Univ.) 27 p CSCL 20K 


G 3/39 


N92-13463 


Unci as 
0053401 


Carlos A. Felippa and Carmelo Militello 
University of Colorado 
Boulder, Colorado 


November 1991 


Prepared for 

Lewis Research Center 

Under Grant NAG3-934 


NASA 

National Aeronautics and 
Space Administration 




developments in variational methods for 
high performance plate and shell elements 


CARLOS A. FELIPPA 
CARMELO MILITELLO 


Department of Aerospace Engineering Sciences 
and Center for Space Structures and Controls 
University of Colorado 
Boulder, Colorado 80309-0439, USA 


abstract 

Hint performance element, are ain.pl. Unit, element, contracted to d .li.t ..^ 
Edition a variant of the ANS tormnlation, died tke 

ssssrci “ - 

illustrate the use of the new formulation. 

1. INTRODUCTION 

For 25 years researcher, have tried to construct ‘best' finite models^ 

problems in structural mechanics. The qu«» !Hh« elements did 

when high order displacement elements dominated the headlines. 

Tot demote the marketplace. The overwhelming preference of finite ekm«t code «ers 

has been for simple elements that deliver engineering accuracy tut coa 

will be collectively called high performance elements, or HP elements for short. 

1.1 Attributes of HP Elements 

Accroaching that general goal gives rise to a myriad of more concrete requirements, 
which are supposed to be addressed in higher or lesser degree during element development. 

Suel p« “-r 

a enneeqaenc. of trym, to acli.v. eatiefacton- accncacy ed. nklmf meal... Bn. o.ler 
items listed in Table 1 call for some explanation. 
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Table 1- Target Requirements for High-Performance Elements 


• Simple: few freedoms, all physical, preferably at comers only 

• Convergent 

• Frame invariant 

• No locking 

• Rank sufficient: no spurious modes 

• Balanced stiffness: not too rigid, not too flexible 

• Stresses as accurate as displacements 

• Low distortion sensitivity 

• Mixable with other elements 

• Economical to form 

• Easily extendible to nonlinear and dynamic analysis 

• Effective local error estimator for mesh adaptation 


The first and foremost requirement is that the element be as simple as possible. 
This is in sharp contrast to the ‘baroque FE period* of 1965-1975 that landed luxuriantly 
ornate elements and culminated with the development of very complex models, including 
elements with nonphysical degrees of freedom. One source of this retrenchment ha. been 
feedback from users of general-purpose finite element programs. As use of these prograns 
expands to more engineers without deep knowledge of “what’s inside the black box* the 
trend in finite element model construction ha. veered toward, the ‘simplest elements that 
will do the job* Further impetus is provided by the gradual realisation that high accuracy 
of complex elements in linear elastoetatics does not necessarily carry over to dynamic and 

nonlinear analysis. , , , 

The balanced stiffness requirement also deserves comment. It follows from the goal 
of attaining reasonable accuracy with coarse meshes. This is illustrated in Fig. 1, which 
shows a convergence study of a classical model problem: the bending of a simply-supported 
square plate under a concentrated central bad. The mesh contains 2 x N x N triangles 
over a plate quadrant. A target ‘accuracy band* of ±1% is taken, somewhat arbitrarily, 
as representative of engineering accuracy for this rather simple problem. The convergence 
characteristics of several triangular elements are taken from the extensive study reported 
in Ref. 2. Although most elements converge, some are too stiff while others an » tool lexib , 
and generally do not enter the accuracy band until the mesh is fairly refined (N _ 8). On 
the other hand, the results labeled ‘FF\ obtained with a plate element based on the free 
formulation (FF) discussed later, lie within the band for all meshes. 

The balanced-stiffness requirement should not be con fused with fast asymptotic 
convergence for fine meshes. Simple elemen ts cannot effectively compete with higher or- 
der elements in this regard, and are not effective in applications that demand very high 
accuracy. What is important is how good art the results for coarse meshes. 
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Figure 1. Convergence study of several plate bending triangular elements 
u reported in Ref. 8. The FF results are from Ref. 8. 


Ti 


L.2 Constructing HP Elements 

The search for high-performance (HP) elements began seriously in the mid 1970s and 
>y now it represents an important area of finite element research in solid and structural 
mechanics. Many ingenious schemes have been tried: reduced and selective integration, 
incompatible modes, mixed and hybrid formulations, stress and strain projections, the free 
formulation (FF), and the assumed natural strain (ANS) formulation. Many researchers 
u-e presenly working to develop such elements. The common theme of the investigations is 


Abandon the conventional displacement formulation 


Several techniques used by researchers in their quest to build better elements are 
itemised in Table 2. It may be noted that many of these were introduced over 20 years ago. 
But it is only recently that a concerted effort is made to combine several tools to produce 
HP elements. For example, the present work draws on items 1, 2, 3, 8, 10, 11 and 12 of 
Table 2. 
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Table 2. Tool* of the Tirade 



Technique 

Fear introduced 

1 . 

Incompatible shape functions 

early 1960s 

2. 

Patch test 

1965 

3. 

Mixed and hybrid variational principles 

1965 

4. 

Projectors 

1967 

5. 

Selective reduced integration 

1969 

6. 

Uniform reduced integration 

1970 

7. 

Assumed strains 

1970 

8 . 

Energy balancing 

1974 

9. 

Directional integration 

1978 

10. 

Limit differential equations 

1982 

11. 

Free formulation 

1984 

12. 

Assumed natural strains 

1984 


1.3 Objective of Present Work 

This paper is part of a series (Refs. 9-12, 15-16) that studies how several HP ele ^ 
construction method, can be embedded within an extended variational framework that uses 
parametrised hybrid functionals. Particular attention is focused on merging the 

items m lan of attack for this unification is flowcharted in Fig. 2. Box connec- 

tion. indicated with dashed lines are not dealt with in the present work. The variational 
extensions, shown on the left of Fig. 2, involve parametria ation of the conventional elastic- 
ity functional, and treatment of element interfaces through generalisations of the hybrid 

aPPr0a ?hf effwtivT"c 2 o 3 nstruction of HP element. reMes on device. 

called •tricks' or Variational crimes that do not fit a prion in the classical variational 
fr^tworic The range of trick, spans innocuous coflocation and finite difference constraints 
fo m^etastic rentes such as selective integration. Despite their unconventional na tore 
tricks are an essential part of the construction of high-performance elements. Collecti ely 
they* represent*!* fun-and-games ingredient that keep, the derivation of HP finite elements 

4 *^The g pre«nt y ?reatment -decriminalises' kinematic constraint trick* by adjoining 
Lagrange multipliers, hence setting out the ensemble on proper variational foundations. 
Placing^ formulations within a variational framework ha. the great ^vantage 
the general structure of the stiffness matrices and foremg vectors of high performance el- 
uents, and of providing theoretical coherence for the systematic derivation of classes 
elements by a combination of techniques. 
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Figure 2 Program of attack for variational formulation of HP element# 


2. THE ELASTICITY PROBLEM 

Consider a linearly elastic body under static loading that occupies the volume V The 
body is bounded by the surface S, which is decomposed into S : S£>S t . Displacements are 
prescribed on S d whereas surface tractions are prescribed on S t . The outward unit normal 

on S is denoted by n s fig. , - _ 

Tie three unknown volume fields are displacements u = tig, in ni 8 # 

e = to, and stresses <r = *y. The problem data include: the body force field b = b< m V, 
prescribed displacements d a l on S d , and prescribed surface tractions t-U on S«. 

The relations between the volume fields are the strain-displacement equations 


e = |(Vu + V r u) = Du 


or 


*a=\ Ki + **/.•) ‘ mV > 


(i) 


the constitutive equations 

^ = Ee or <Jij = Eijki*ki ** 

which will be assumed to be invertible, and the equilibrium (balance) equations 
— div <r = D'ff = b or o'.-y.y + 6. = 0 “ 


( 2 ) 

(3) 
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in which D* = -div denotes the adjoint operator of D = §(V + V 7 ). 

The stress vector with respect to a direction defined by the nmt vector v is denoted 
as = «r.v, or <r w = <r i} Vj. On 5 the surface-traction stress vector is defined as 

9 n - 9 XL, Or <Tni = <*</»/. ( 4 ) 

With this definition the traction boundary conditions may be stated as 

* £ or «*/»/ = & on 5*, ( 5 ) 

and the displacement boundary conditions as 

U » 4 or Ui ■ i on S+ ( 6 ) 


5. NOTATION 


3.1 Field Dependency 

In variational methods of approximation we do not work of course with the 
field, that satisfy the governing Eqs. 1-3 and M, but with tndependent Mg; 

which are subject to variations, and dependent (secondary, associatedderived fields, wh*h 
are not. The approximation is determined by taking variations with respect to the mde- 

^ ndCn kntring the notation introduced in Refs. 9 and 10 an iWependent/yeamdfield 
will be identified by a superposed tilde, for example <L A dependent field is identified by 
writing the independent field symbol as superscript. For example, if the displacements are 
independently varied, the derived strain and stress fields are 

e « = 1(V + V T )u = Dii, = E«“ = EDu. (7) 

An advantage of this convention is that n, e and » may be reserved for the exact fields. 


3.2 Integral Abbreviations 

Volume and surface integrals will be abbreviated by placing domain-subscripted 
parentheses and square brackets, respectively, around the integrand. For example: 

(/ )v Wjf/jv, |/|,« f s t*s. VU.-'fjds. [/I, IjdS- m 

If f and g are vector functions, and p and q tensor functions, their inner product over V u 
denoted in the usual manner 

(f, g) v 12 J f.g dV = jf fi 9 i dV, (p, q)v =' ^p.qdV = ^ <&• ( 9 ) 

and similarly for surface integrals, in which case square brackets are used. 

3.3 Domain Assertions 


The notation 

(a = b)vt [ a = ^5» [ a ** I a “ 

is used to assert that the relation a = b is valid at each point of V, S, S d and S t , respectively. 
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Figure 3. Internal Interface example. 


3.4 Tntt™* 1 Tnterfacca 

In 54-5 we construct hybrid ea national principles in which boundary 
m .< ££ZS* I * iS*. of boundary to tactad. 

daaignatad aa S|. Tbu* S;J , uSuSl . («) 

Oa 5, nai.bav diapUa.mani. .or taa.ta. « P"»f>rf. V""* 1 ! 

“i-.d 

^ i^/M .n.oo.h point, of St aba oni. nonnal. »♦ and tf pout. » 

OPPO ‘ , ^“' abbreviation. ot Ev- M ganaraliaa aa follow., o.iag Fig- 3 for 
nesa. A volume integral ie the aum of integral* over the eubvolumea: 

(/)v « (12) 

An integral over 5* include* two contribution*: 

|,| s .tf/ s s *dS + / sr p-d5. 

.h, r « ,+ and daaotaa tha vain, cf tb. inland , on St and ST. 

„0 valna. may ba difarant if , i. dia.ontinuon. oa involva. a pro,actto» »» *• 

The appearance of S v is a natural coneequence of use of finite element* with duoon- 
tinuou* displacement*. FoUowing a finite element discretisation, the union of interelemen 


(13) 
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boundaries becomes S { . This boundary is generally nonphysical because it depends on the 
discretisation. 1 * * 


4. THE ELASTICITY FUNCTIONALS 

The variational principles of linear elasticity are based on functionals of the form 

n - V - P, ( 14 ) 


where U characterises the internal energy stored in the body volume and P includes other 
contributions such as work of applied loads and energy stored on internal interfaces. We 
shall call U the generalized strain energy and P the forcing potential. 

The functionals considered in this section include independently varied displace- 
ments. The class of dual functionals such as the complementary energy are briefly cover'd 
in §5.5 for completeness, but are not required in the finite element developments of §6ff. 

4.1 Generalised Strain Energy 

The generalised strain energy has the following structure: 


U = \jn{*i**)v + J12 (*■ *)v + J13 (*. «“)v + » e V 

where ji\ through 3 $ 3 are numerical coefficients. For example, the Hu-Washisu principle 
is obtained by setting J 13 = -1, Jis - h ha = I, other * bein * ,ero * The^ matrix 
representation of the general functional Eq. 15 and the relations that must exist between 

the coefficients are studied in §5.1. 


4.2 Hybrid Forcing Potentials 

Variational principles of linear elasticity are constructed by combining the volume 
integral of Eq. 15 with the farcing potential P. Two forms of the forcing potential, called 
pd ga j p* ja the sequel, are of interest in the hybrid treatment of interface discontinuities. 
The ^-generalised (displacement-generalised) forcing potential introduces, as anticipated 
in §3.4, an independent boundary displacement field d over 5<: 

P d { u, i, d) = (b,u)v + [J„,u - 3] 5 , + Ms. + [2r„,u - d] Sj . (16) 


The t-generalised (traction generalised) forcing potential introduces an indepen- 
dently varied traction displacement field t over 

P*(u, J,t) = (b,u)v + [t,u- d]s 4 + Mis. + Mk- (I 7 ) 


The ‘conventional* form P° ot the forcing potential is obtained if the interface 
integral vanishes and one sets [t = <r„] s . If so f* and P d coalesce into P e , which retains 
only two independent fields: 

P*(u, i) = (b,u)v + - d] 5< + MI*.* (18) 


1 If there are physical internal interfaces — for example a sudden thickness or material change — 

it is common practice to select the mesh so that these natural interfaces are also interelement 

boundaries. 
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4.3 Modified Forcing Potentials 


Through various manipulation* and assumptions detailed in Ref. 10 the forcmg 
potential P 4 * may be transformed to 

P^u, i, d) = (b,u)v + Ms, + [#„,u - d|s. ( 19 ) 

where the all-important surface dislocation integral is taken over 5 rather than Sj. One of 
the assumptions is that displacement boundary conditions, Eq. «, are exactly satisfied on 
S d . This expression of P* is used in the sequel A similar technique can be used to adjust 
P‘, but that modified formula will not be required in what follows. 

4.4 Complete Functionals 

Complete elasticity functionals are obtained by combining the generalised strain 
energy with one of the forcing potentials. For example, the d and t generalised versions of 
the Hu-Washiiu functional are ______ 

n£r » Uw - P*, T§y**Uw-P*' t 20 ) 

where U w is obtained by setting jn =» ;'is * L “ -1 > othtn sero, m Eq. 15. 


5. MATRIX REPRESENTATION OF ELASTICITY FUNCTIONALS 
The generalised strain energy of Eq. 15 can be presented in the matrix form 


U 


i/ v (* ^ 


in 

symm 


in 

722 



( 21 ) 


The symmetric matrix 3 



in 

in 

in 

J = 

ill 

hi 

hi 


w ii s 

in 

in. 


( 22 ) 


characterises the volume portion of the variational principle. Using the relations = Ee, 
= EDu, e' = E -1 <r, and e" = Du, the above integral may be rewritten m terms of the 

independent fields as 




(J e u) 


jnE 1 
Jial 
JiaD r 


JijI 

i . wE 

jnsD r E 


j laD 
?D 

jisD r ED 



dV. 


(23) 


a To justify the symmetry of J note, for example, that /is(*,«*)v - $/ia(<r,«“)v + 

Jiis(*',» m )v.»adsoon. 
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5.1 


First Variation of Generalised Strain Energy 
The first variation of Eq. 15 may be presented as 

SU = (Ae, Sc)v + (A*, 6*)v ~ (<Hv J . + Wn> * u U» 


where 


Ae = jne' + j\ 2 « + Jis®“» 
Ae *» jut + jn** + Jase", 
■ Jis* + h**" + *»*"• 


( 25 ) 


The last two terms combine with contributions from the forcing potential variation. 
For example, if P = P* the complete variation of II e * U - P* is 

5TI* = (Ae, 8i) v + (Ae, Si) v - (dW + b, 6u)v + K - - [5 - d, »«] j 4 - (26) 

Using P* or J* does not change the volume terms. The first variations of IT* and 
II* are studied in Refs. 9-11 for a more restrictive class of functionals, namely IL,. Th 
Euler equations associated with the volume terms 

Ae>0, Ae = 0, dlvV + b = 0, (27) 

are independent of the forcing potential A “weighted residual* inte^retationof^27in 
terms of the field equations is given in §5.4. For the moment we notethat forconsistency 
of the Euler equations with the field equations of §2 we must have Ae * 0, Ae - 0 and 
J = 9 if the assumed stress and strain fields reduce to the exact ones. Consequently 


iii + ii2 + in = o» 
iia + ha + J» “ °» 
ilS + J29 + m 1* 


(28) 


Because of these constraints, the maximum number of independent parameters that 
define the entriee of J U three. 

5.2 Specific Functionals 

Expressions of J for some classical and parametrised variational principles of elas- 
ticity tabulated below. The subscript of J is used the identify the functionals, which 
are Led roughly in order of ascending complexity. The fields included in parentheses after 
the functional name are those subject to independent variations m . 

Potential energy (u): 

0 0 0 
0 0 0 
0 0 1 


(29) 


Stress-displacement Reissner, also called Hellinger-Reissner, (8, ft): 




-10 1 

0 0 0 
10 0 


(30) 


Unnamed stress-displacement functional listed on p. 116 of Ref. 18 (*,u): 


Jt; = 


10-1 

0 0 0 

-10 2 


(31) 
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Figure 4 Graphical representation of the Jafii functional* 
Strain-displacement Reissner-type as listed on p. 116 of Ref. 18 («,u): 


Hu-Washisu 3 (5,e,u): 


'0 0 0 

J s = 0 -1 1 - 

0 10 


’ 0 -1 1 
J w — “1 10. 

10 0 


One-parameter stress-displacement family (>,u) that includes Up, Ur and Uu as 
special cases (Refs. 8-10) 

— *# 0 nr 

J, = 0 0 o . ( 34 ) 

[l 0 1-Tf. 

One-parameter strain-displacement family (e,u) that includes U P and U s as special 
cases (Ref. 9) q 

J, = 0 —0 P ( 3S ) 

0 p l-p. 

» There are several functional, that carry this name, transformable from tone! to th ” ugh 

integration by parts. That corresponding to Jw » the third form listed in §2.3 of Ref. 24. 
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Two-parameter family (»,i, *) that include* U, and Cf 7 a* special case* (Ref. 9): 

J* - (1 - flJ-r + (1 “ 1 ) 3 * -{l-P-Wr 

’-7(1-0) 0 7(1-0) 1 (36) 

0 -0(1-7) 0(1-7) 

7(1-0) 0U-*7) 1-0-7 + 20'lJ 

Three-parameter (a, 0, 7 ) family (*, *,&) that include* U w *nd U fl a* «pecial ca*e« 
(Ref. 9): 

J ttfli “ aJnr + (1 - _ 

-7(l-0)(l-«) ~ a a + 7(l-0)(l- a ) (37) 

( -a a-0(l-7)(l-a) ? 1 “*M l -|S ) , * 

i a + 7( 1 -0)(l-«) 0(1-7)(1-®) (l-0-7 + *07)(l-«)J 

The last form, which contain* three independent parameter*, *uppUe* all matri- 
ces J that satisfy the constraint* of Eq. 28. It yields stres^pl^ement functionaJs for 
a a 0 s 0, strain-displacement functional* for o = 7 - 0. “J three-field (»tie*»-s - 

displacement) functional* otherwise. A graphic representation of tlu* functional m (a, 0, 7 ) 
apace is given in Fig. 4. 

The specialisation 7 = 1 , 0 - 0 of 3 a01 is of particular interest: 


a — 1 -a 1 

—a a 0 

1 0 0 


(38) 


The associated functional IT a might be called the -gencralued- Hu- Waslusu fun^ 
tional since it reduce, to 3 W for a * 1. But because of it. special rel^mnwith the ANDES 
formulation covered in §8-11, n a will be herein referred to as the ANDES functional. 

5.3 Energy Balancing 

A prime motivation for introducing the 3 coefficients as free parameters is opti- 
misation of finite element performance in the balanced-stiffness sense of Table 1. T e 
determination of “best* parameters for specific elements relies on the concept of energy 
balance. Let U(e) = £(Ee,<)v denote the strain energy associated with the steam field t. 
If E is positive definite, I/(e) is nonnegative. We may decompose the generalised strain 
energy into the following sum of strain energies: 

U = u( e «) + Wl U(e* - e) + ^{e - «“) + «flsl/(e" - «'), (39) 

where Up(« u ) - Up is the usual strain energy, and 4 
tDl 3 J(/u+j32— J33 + 1), U»3 = |(-Jll+il3+ft3-l)« <"S l“ \ [jll ~072 + J33~ l)« (*°) 


Eq. (39) is equivalent to decomposing J into the sum of four rank-one matrices: 


J = 


0 0 0 
0 0 0 
0 0 1 


+ U»1 


1 -1 0 

-1 10 

0 0 0 


+ u>3 


0 0 0 
0 1-1 

0-1 1 


+ u>3 


1 0 -1 

0 0 0 

-10 1 


(41) 


Decompositions of this nature can be used to derive energy-balanced finite elements 
by considering element -patches* under simple load systems. This technique is discussed 
for the one-parameter functionals generated by J, in Refs. 6, 8-11. It is impor an o no 
that the j coefficients may vary from element to element. 

4 A* shown in §5.4, these coefficients may be Interpreted as field-equation-residual weights, 
hence the notation. It is conjectured that for stability the j coefficient, should be confined 
80 that Wi > 0, but this remains to be proven. 
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5.4 


( 42 ) 


Interpretation of Euler Equations 
Eq». 27 gain physical meaning if they are rewritten u 

A« » wi(e' - i) + «*(«* - «') = °. 

A# * v>i{* — **) + — # *) “ 

dhr s' « dhr [*“ + ^{v* - c“) + «*(*• - *)] * -b, 

, M D. t E“ l * - 4 = 0 as well a* i - •* ■ 

-r*- *•■ v®tr rs£: 

iuiSi *. Mi’i— a. - o . -a,' 

kinematic and constitutive field equation*. On the other hand, dlv + 

combination of the equilibrium equation*, Eqs. 3, and the other . 

If the i coefficienU are auch that a weight vanishes (aee abo footnote 4}, that partK 

^ . , . (■ thsEnler equation* and must be viewed aa being satisfied 

ular field equation drop* out from the EuJer equation* ana u» „ , j the 

Points the way for constructing U of Eq. 15 by the method of weighted residuala. 

5.5 Functionals without Independent D isplacements 

a.** “S 

to displacement-derived strains: 

Uu = faf + jnS + lisa'W - (dhr^.njv -K.«U- (43) 

If we now assume that the equilibrium equation, div #+b = 0 and traction boundary 
conditionT ~ tMft a priori,^ may be dropped and w. are left with the ,«»eruW 
complementary energy functional 

U — U c -= fj'u(e, *?)v + Jia(*. «)v + a 8 V ■ ^ 

Taking account of the a priori conditions, the first variation becomes 

6U e = On®* + J 1 3 ® + ®*» 6 *) v + (i™* + Se ^ v ’ 

and for consistency we must have j\i + jia = -1. Jia + ha = j* f°il° wa 
represented as inThe matrix form of Eq. 21 with a J that depend, on a single parameter 


J,= 



(46) 


Here p — 0 gives the classical principle of total complementary energy whereas p - 1 
gives the functional N(i,e) listed on p. 117 of Ref. 18. 

0. FINITE ELEMENT DISCRETIZATION 

In this section assumptions invoked in the finite element discretisation of the fu “ c ‘ 
tional IT* for arbitrary J are stated. Following usual practice in finite element work, the 
componentsof street and strains are arranged as one-dimensional *ray. 
moduli in E are arranged a a a square symmetric matrix. In the seque , 

noted, we consider an individual element of volume V and surface 5 : 5 t U S d U S„ where 
Si is the portion of the boundary in common with other elements. 
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6.1 Boundary Displacement Assumption 

The boundary displacement assumption is 

[d * N*v]s. 


( 47 ) 


Here matrix N„ collects boundary shape functions for the boundary displacements 
d while vector v collect, the “visible- degree, of freedom of the element, alw called ^the 
connector.. These displacement* must be unique on common element boundaries. Tbu 
continuity condition is met if the displacement of a common l>°uu<l«y portion is wuvdy 
specified by degrees of freedom located on that boundary. There are no derived fields 

associated with <L 

6.2 Internal Displacement Assumption 

The displacement assumption in the interior of the element is 

(u = N.q)vr, ( 48 ) 

where matrix N. collect, the internal displacement shape functions and vector q collects 
generalised coordinate, for the internal displacements. The assumed u need not be contin- 
uous across interelement boundaries. The displacement derived fields are 

(«• = DNq * Bq)v , (#* = EBq)v • ( 49 J 

To link up with the FF and ANS formulations, we proceed to break up the internal 
displacement field as follows. The assumed u is decomposed into rigid body, constant 
strain, and higher order displacements: 

u = N r q, + N.q. + N^q*. ( 50 ) 

Applying the strain operator D = |(V + V r ) to u we get the associated strain field: 

e“ = DN r q, + DN^ + DN^q* = B r q, + B c q« + B^q*. (51) 

But B r = DN r vanishes because N, contains only rigid-body modes. We are also 
free to select B 0 = DN C to be the identity matrix I if the generalised coordinates q e are 
identified with the mean (volume-averaged) strain values e*. Consequently Eq. 51 simplifies 

l ° e -«?‘ + eX-e , ‘+Bhq k , ( 52 ) 


in wmen , . _ f*o\ 

q, 5l* = (e“)v/v> ( B *)v = 0. ( 53 ) 

where v = (l)v is the element volume measure. The second relation is obtained by integrat- 
ing both sides of Eq. 52 over V and noting that q* is arbitrary. It says that the mean value 
of the higher-order displacement-derived strains (also called the devtatonc displacement- 
derived strains) is »ero over the element. 

6.3 Stress Assumption 

The stress field will be assumed to be constant over the element: 

{i=e)v- ( 54 ) 

This assumption is sufficient to construct high-performance elements based on the 
free formulation (FF) developed in Refs. 4-5. As discussed in Ref. 11, the inclusion of 
higher order stress modes (deviatoric stresses) iin Eq. 56 is computationally effective if 
these modes are divergence free, but such a requirement makes extension to geometrically 
nonlinear problems difficult. The only derived field is 

(e* =r = E -1 tf)v. -( 55 ) 
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6.4 Strain Assumptions 

The assumed strain field i is split into a mean constant strain 7 and a higher order 

variation (the devi&toric strains) i 

(e = ?+e,j = e + Aaa)v, ( 56 ) 

where 7 - (e) v /«, matrix A* collects deviatoric strain modes with mean sero value over 

the element: /r7A 

(A*)v — 0, ( 57 ) 

and a collects the corresponding strain mode amplitudes. The only derived field is 

(*• = Ee = Ee + EA^aJv . ( 58 ) 


T. UNCONSTRAINED FINITE ELEMENT EQUATIONS 

For simplicity we shall assume that all elastic moduli in E are constant over the 
element. Inserting the above assumptions into IT 1 with the modified forcing potential 
of Eq. 19 we obtain a quadratic algebraic form, which is block-sparse on account of the 
conditions stated in Eqs. 53 and 57. Rendering this form stationary yields the finite element 
equations 


JuvE- 1 

j l2 vl 

0 

-p? 

y l3 vi 


0 

0 

0 

0 


0 

-Pr 

0 

0 

0 

iisvl - P« 

jisul 

0 

0 

-Pa 

0 

J13H* 

0 

L 

0 

0 

0 



where 

K*, = (AJEA<)v - R = (B^EAa)v, 

L = (NL]s, Pr - INfjs, P« = (NLls, Pa = lNL]s, («°) 

f r = (NfbV, C, = (N?b)v, fa = (N^bjv. U = [Njt] s „ 

in which denotes the projection of shape functions Nj on the exterior normal n, and 
similarly for N r , N e and N*. Those coefficient matrix entries that do not depend on the 
j coefficients come from the last boundary term in Eq. 19. 


7.1 The P Matrices 

Application of the divergence theorem to the work of the mean stress on e u yields 

(7, e“)v = (7, 5“ + B„qjv - v* T ?‘ + * T (Ba)vq* = u7 r «* (61 ) 

= [^n,«] s = [^.Nrq, + N e r +N fc q*] s = ^(Prq, + Pci“ + PaQa)- 
Hence P r = 0, P c = vl, P h = 0, and the element equations simplify to 


juv'E 1 

j x2 vl 

0 

0 

(J13 - 1)1 

iiawl 

jnVE 

0 

0 


0 

0 


0 

0 

0 

0 

0 

0 

0 

(j ’13 — 

323 vl 

0 

0 

3$ 3^E 

0 

0 

32 3®* 

0 

0 

L 

0 

0 

0 

0 


0 

l t i 


7 ' 


r 0 

0 

0 


e 


0 


0 


a 


0 

0 

0 


Or 

1 = < 

Ur 

0 

0 


e* 


fgu 


0 


q h 


Uk 

0 

0 . 


a V , 


< f V 4 


( 62 ) 
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The simplicity of the P matrices is essentially due to the me an- plus-de viator splitting 
of Eq. 52 for e“. If this decomposition is not enforced, P r = 0 but P, - [B e ) v - vB 0 

Pfc ■ (Bfc)v* 


8. KINEMATIC CONSTRAINTS 

The ‘tricks’ we shall consider here are kinematic constraint* that play a key role m 
the development of high-performance FF and ANS elements. The* are matrix relation, b^. 
tween kinematic quantities that are established independently of the variational equations. 
Two types of relations will be studied. 

Constraints Between Internal and Bounda ry Displacements 

Relations linking the generalised coordinate, q of Eq. 48 and the conn^tors v were 
introduced by Bergan and coworkers in conjunction with the free formulation (FF) of fim 
element, (see Ref. 5). For simplicity we shall assume that the number of freedoms ^^and 
q is the same; removal of this restriction is studied in Ref. 11. By collocation ofuat the 
element node point* one easily establishes the relation 


8.1 


v = Grq, + G c q* + G*q* * Gq, 


( 63 ) 


where G is a square transformation matrix that win be assumed to be nonsmgular. On 
inverting this relation we obtain 


i-i 


Hv, 


or 



(64) 


The following relations between L (defined in Eq. 80) and the above submatrices 
hold as a consequence of the individual element test described in |9.3: 


L r G r * 0, L t G, = vl, uH« 


(65) 


If the splitting of Eq. 52 is not enforced, however, the last two become 

L T G e = vB e , P e H* + PfcHa =» L t or P e = L T G e . (66) 

Since P. = oB 0 , these relations coalesce (see Ref. 5). 

8.2 Constraints Between Assumed Deviatoric Strains and B oundary Displacements 

Constraints linking e to v are of fundamental importance in the assumed natural 
strain (ANS) formulation. The effect of these constraints in a variational framework is 
analysed in Refs. 15 and 16. In the present study we depart from previous work mjhat 
only the deviatoric strain t, e Jt are assumed linked to v mea " *J ra ™ e are 

obtained variationally. Consequently we shall postulate the foDowmg relation between 
assumed deviatoric strain amplitudes and nodal displacement connectors: 

a = Qv, (67) 


where Q is generally a rectangular matrix determined by collocation, least squares or 
other fitting methods. An example showing the construction of Q is given m §11.4. The 
individual element test described in §9.3 requires that Q be orthogonal to G r and G e . 

QG r = 0, QG„ = 0. (68) 
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8.3 Limitation Principle* 

Strain sumption* made concurrently with displacement «,umptioM areconfin^i 
, .. . . ta4 . -nnciole* similar to those stated by Ftaeijs de Veubeke for stress-displacement 

i^ixed elements (Ref. 13). This issue was discussed in Ref. 15 for a more restricted strain- 
r p £,»“. “ Ura formulation. UmiuUn. pri~ipU. for tlu formula, pr- 

tented here remain to be studied. 

0 . VISIBLE STIFFNESS EQUATIONS 

On enforcing the constraints a = Qv, q,. * H r v, q* == H«v * v J* v * 
q^ = H fc v, through Lagrange multiplier vectors A*, A r , A c , and A*, respective y, we ge 

the augmented finite element equations 


AivE' 1 Ji7vl 0 

Jiavl J»vE 0 

0 0 o 
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0 0 
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0 Ha 


q t h T tr l i . Hi 0 
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5* 




Qk 

► =r ■ 



A* 
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Ar 
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A. 


0 


Xk 


0 


L V , 




(69) 


Condensation of all degrees of freedom except v yields the visible 8 element stiffness 
CqUati ° M Kv = (K* + Kj»)v = f, (70) 

Wh< " K* = v" 1 LEL r , ( 71 ) 

K k =* JJ3 hZK'hBh + + Q T H, r Hfc) + *aQ r K a 4Q, (72) 

f = f;+Hff,r + u-‘i% +nr*t. ( 73 ) 

Following the nomenclature of the free formulation, we shall call K* the basic stiff- 
ness matrix and Ka the higher order stiffness matrix. 

9.1 Relation to P revious HP Element Forffiqhtfeafi 

If J = J, of Eq. 33, - 1 - 7, 372 = *» - , *»d we recover the scaled free 

formulation stiffness equations considered in Refs. 6, 8 and 10: 

K* = (1 - 1) 1-1>0. (74) 

On the other hand, If we take J = J a as given in Eq. 38, jw - <*. Jos = - 0 and 

"° bUin K fc a aQ r KoaQ. «>°, < 75 > 

which is similar to the stiffness produced by the ANS hybrid variational 

in Refs. 15-16, in which the forcing potential P* was used instead of P d . The variant of ANS 

» The qualifier visible emphasise, that these are the stiffness equations other elements W, 
«d consequently are the only one. that matter in*>far as computer implementation on a 

displacement-based finite element program. 
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considered herein will be called the assumed natural deviatoric strain (ANDES) formulation 
in the sequel. The name is apt in the sense that what is being assumed are deviatonc rather 
than total strains, and that this assumption only affects the higher order stiffness. 

But the term with coefficient in Eq. 72 is new. It may be viewed as coupling the 
FF and ANDES formulations. It is not known at this time whether Eqs. 70-73 represent 
the most general structure of the visible stiffness equations of HP elements. 

9.2 Recovery of Element Fields 

For simplicity suppose that the body forces vanish and so do f, r , and °® 
account of Eqs. 60. If ▼ is known following a finite element solution of the assembled 
system, solving Eqs. 69 for the internal degrees of freedom yields 

«= v- x L t v, 7 = Ee, a= Qv, q, = H r v, e“*7, q* = H h v, ^ 
A.-(»I^Q + * 3 R r Ha)v, Xr = 0, A.-0, A fc = (*sHQ + *sK ffc H„)v. 

It is seen that the mean strains 7, e* and e* = B -1 7 coincide, and of course so do 
the mean stresses. But if the body forces do not vanish the mean stresses and mean strains 
recovered from different fields will not generally agree. 

It is also worthwhile to note that a nonsero Lagrange multiplier vector flags a devi- 
ation of the associated fields from the variationally consistent fields that would result on 
using the unconstrained Eqs. 62 without "tricks* . 

9.3 The Individ" *1 Element Test 

To conclude the general formulation we investigate the conditions under which HP 
elements based on the foregoing seeing pass the individual element test of Bergan and 
Hanssen described in Refs. 3-6. To carry out the test, assume that the free floating” 
element* under sero body forces is in a constant stress state *o, which of course is also the 
mean stress. Insert the following data in the left-hand side vector of Eq. 69: 

7=a 0 =y, e = E -1 Oo. ** = °, q, - arbitrary, e“ = e“ = E“ l 7 0 , q* = 0, 

A„ = 0, A r — 0, A„ = 0, A* = 0, v=»G r q r -l-G«e ,, = G r q r + G e E- 1 »o. 

Premultiply by the coefficient matrix, and demand that all terms on the right-hand 
side vanish but for f, = L<r 0 . Then the orthogonality conditions listed in Eqs. 65 and 
68 emerge. This form of the patch test is very strong, and it may well be that relaxing 
circumstances can be found for specific problems such as shells. 

10. DISCUSSION 

At this point it is useful to recapitulate key points of the previous development, wd 
to connect this material with some of the techniques of Table 2. The chief property of HP 
elements constructed with present methods is the decomposition of the element stiffness 
equations displayed in Eq. 70; a property that of course subsists at the assembly level. 

The basic stiffness matrix has a universal character as no coefficients j appear m 
Eq. 71, clearly K b is independent of specific variational principles. Given the constant stress 
state introduced in Eq. 54, K* depends only on the assumed boundary motions. It can be 
constructed (and programmed) once and for all for each element type. As emphasised m 
Ref. 5, the main function of Kb is to provide convergence. 

The higher order stiffness given in Eq. 72 serves two other functions: stability and 
accuracy. The basic stiffness is generally rank-deficient 7 because its rank cannot exceed 


* Mathematically, the entire element boundary is tractlon-specifled, S 3 S„ 

1 Except In simplex element®, for which K = K&. 
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that of E; thus a key function of K* is to stabilise K by raising its rank to the correct one. 
The second function, which has gained importance in recent work, is to increase solution 
accuracy for coarse grids. Here is where the j coefficients play the important role noted 
in §5.3. These coefficients may vary from element to element, despite the fact that this 
variation implies that the variational principle changes from one element to another. Thus 
the •element mixability’ requirement of Table 1 is fulfilled without tears. 

10.1 The Free Formulation 

The present methodology was initially pursued to justify variational* the original 
FF h = 0) of Ref. 5 as well as the scaled FF (7 * 0) of Refs. 6-8. Thus it is not 
surprising that those element construction techniques fit naturally in the present variation 

faLJt by .Imply UUM J = I.- “ * “T * 

remove the restriction that the dimension of vectors q and ▼ be the same. One of 
techniques advocated to allow dim(q) > dim(v) involves extending Eq. ! 54 with 1 deviatonc 
stress assumptions, and thus requires a generalisation of Eqs. 59 and 62. Whether sn 
generalisation is practically worthwhile is unclear at this time. 

10.2 The AN3 Fo rmulation 

The conventional ANS formulation as presented in Refs. 1 and 19 constructs total 
strain fields e (not necessarily integrable into displacements u“) gaged through generalu 
strain coordinates a as e « An. These coordinates are eventually linked to the connectors 
v via matrix expressions of the form a = Qv, leading to an element stiffness of the form 
K = Q r K-Q, where K* is the generalised stiffness in terms of *. The restriction to 
deviatoric strains in § 6.4 is motivated by two interrelated factors: (a) the stram-assumed 
stiffness "flows’ to the higher order stiffness, where it can be naturally scried by using 
J = J„, and even intermixed with FF contributions as Eq. 72 shows; and (b) the basic 
stiffness’ of the element, derived separately, can be used to insure convergence. 

10.3 Projectors and S/R Integration 

The so-called *B-bar* approach is based on expressing the element strains as* 

e = Sv (78) 

where B, which cuts off the “harmful’ portion of B*. is constructed by various ad-hoc 
devices such as strain projection, selective and/or uniform reduced integration. These 
time-honored schemes are well covered in Ref. 14. They are easily included in the present 

setting if B admits the decomposition 

S = 5 + AaQ, ( 7Q ) 

where Q is not position dependent and I - Ev provides the mean strains, which are 
discarded in favor of Eq. 76. This decomposition can be usually earned out in several 

way*. 

11. EXAMPLE: A 9-DOF ANDES PLATE BENDING TRIANGLE 

The first element constructed with the ANDES formulation is a three-node Krchboff 
plate-bending flat triangle with the usual nine degrees of freedom. The derivation is briefly 
covered here as it illustrates the essential steps in forming the higher order stiffness of such 
elements. These steps are outlined in “recipe’ form in Table 3. This Table "states the 
arguments of § 6.4 in a more physically oriented sense, which is closely aligned with the 

terminology of Ref. 19. 

• This is a slight variation from the usual notation, necessitated by the use of the single overbar 
to denote average or mean value*. 
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Table 3 


Construction of K* by the ANDES Formulation 


SUp 1. Select reference lines (in 2D elements) or reference planes (in 3D 
elements) where “natural straingage* locations are to be chosen. By appropriate 
interpolation express the element natural strains e in terms of the stramgage 
readings* g at those locations: ^ _ (80) 

where ? is a strain field in natural coordinate, that must include all Want 
strain states. (For bending elemenU the Urn ‘strains* is to be interpreted in a 
generalised sense, viz. curvatures.) 

Step t. Relate the Cartesian strains e to the natural strains: 


e = T? = TA<g = Ag 


(81) 


at each point in the element. (If e = «, or if it is possible to work throughout in 
natural coordinate*, thi* *tep is skipped.) 

sup S. Split the Cartesian strain field into mean (volume-averaged) and devi- 

atone strains: _ /aM 

e = e + e d => (X + Aa)g, l 82 ' 

where X « (TA«W/v, and e d = Adg has mean sero value over V. (This step 
may also be carried out on the natural strains if T is constant, as is the case for 
the element derived here.) 

SUp 4. Relate the natural straingage readings g to the visible degrees of 
freedom g = ^ (83) 

where Q is a straingage-to-node displacement transformation mat^ T^qu« 
by which this is accomplished vary from element to element and is difficult to state 
rules that apply to every situation. In the element derived here Q is constructed 
by direct interpolation over the reference lines. (In general there is no internal 
displacement field u* such that e « Du*, so this step cannot be done by simply 
integrating the field of Eq. 81 over the element and collocating u at the nodes.) 

SUp 5. The higher-order stiffness matrix is given by 

TCft = qQ r K — iQ. where Kad = AjEAj dV, (84) 

where a > 0 is the scaling coefficient supplied by the functional of Eq. 38. 


11.1 n^m.tric Relations 

The triangle has straight sides. Its geometry is completely defined by the location of 
its three corners, which are labeled 1,2,3, moving counterclockwise. The triangle is referred 

the corner coordinates Xi, j* satisfy the relations *i + *a + *s = 0 + . V7 t !*. , * 

Coordinate differences are abbreviated by writing *</ = *< " Xj, etc. The signed tnang e 

area A is given by the formulas 

2A = *3iy31 - *3lV31 = *32yi 2 “ = *13Sto “ *23yi3, (85) 
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- « *+• ' ‘ > °- » * ■* 

&££ £E2££5Xt 82L£ of . otroight-oided trU.,1. » »»Ud= 

ft = ^(xiS/k - *kl/y + *»j«> + ^ 

where t i and Jfe denote positive cyclic permutation, of 1, 2 and 3; for example, t = 2, j = 3 
r-T’Th^fore “ft/3x - y»/iA and 3ft/3y = x k ,/2A. Other intrinsic dnnension. and 

ratios used below are 
the side lengths. 

11.2 Displacements, Rotations, Side Coo rdinates 

Sbc. wo on dealing with o KirekloJ elomeot, ilo diploeomont Seld ia oomP 1 **'!' 

freedom of the element collected in v are 

1*1 0,1 wa 0*3 0»3 W8 ^» 3 1 • 


v r = 


Over the three side. 1-2, 2-3 and 3-1, traveled connterdockwue, we define the 

“ - irsOTss ?*=& 

• _ 1 * —C mkei»wmfa Tnffl 


dx 

dx 

dx 

— .1 ss lie, 
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— — =* X33, 

Op29 

0 Msi 

(89) 

dy 

t = yai» 
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dy 

^- = V33, 
Op23 

dy 

■5 — = yw- 
d/*si 



11.3 ^j^tural Onrvatnrea ... 

be caU^hr^t^^rel 

tion of Eq. 89: 
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(90) 


The inverse of this relation is 
3 2 to 

d? 

93. 


1.x 

[ 4d4 2 




1/33 yi3 
*33*13 
L 1/33*31 + *331/13 


l/silbi 

*31*31 

1/31*13 + *131/31 


S/131/33 

*13*33 



[fit] 

< 

fit 


1ft 


( 91 ) 
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or, in matrix form 


* = Tx- 


- (92) 


11.4 Curvature Sampling 

The reference linet referred to in Table 3 are the three triangle eidee. The natural 
curvature* are assumed to vary linearly over each reference line, an assumption which is 
obviously consistent with cubic beam-like variations of w over the sides. A linear variation 
on each side is determined by two straingage sample points, which we choee to be at the 

corners^ ^ aide the isoparametric coordinates that var X frora ** 

corner t to +1 at corner j. These are related to the Hi, coordinates as = 2 Hi ~ 1 - T “ en 
the natural curvature over side ij is given by the beam formula 


Xij 


d 7 xo 





r ** 1 

IZ 

36,-1 

36, + 1] - 

1 

1 w > 





(93) 


where 9 n denote the rotation about the external normal direction » on side ij. Evaluating 
these relations at the nodes by setting & = ±1 and converting normal rotations to x-y 
rotations, we build the transformation 
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The left hand side is the natural straingage reading vector called g in Table 
so we can express this as the matrix relation 


ti?i 
0*1 
Byl 

0*2 
0 V 2 
W* 

0*3 
0y3 
(94) 

3 and 


Qv. 


(95) 


11,5 Curvature Interpolation 

The six gage readings collected in g provide curvatures along the 3 triangle sides 
directions at two comers. But 9 values are needed to recover the complete curvature field 
over the element. The 3 additional values are th e natur al curvatures at the missing corner. 
We obtain these values by adopting the following rule: Cylindrical bending with linearly 
varying curvature along a side direction is to be exactly represented Another way of stating 
this is that the side curvature x.y » to be constant along lines normal to side ij. This makes 
the element insensitive to bad aspect ratios on -strip bending* if each element has a side 

oriented in the direction of the strip. _ 

To apply this rule consider side 1-2. The natural curvature xia = 3'W*Mi2 
this side is defined at nodes 1 and 2 by the first two rows of Eq. 94. For node 3 take 


Xiala 




Ai3 X12I1 + ^21 X12I3 » 


(90) 
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where A xa and A„ are defined in Eq. 87. Since we now know the value, of Xxa = ia 

at the three corner, we can use the standard linear interpolation over the entire triangle. 

Xia - Xiali ft + Xial, ft + Xuls ft « Xul x (ft + *«ft) + X»la (ft + A «ft)- ( 97 ) 
Proceeding analogously for the other two side, we construct the matrix relation 


Is) 


ft + Aiafc 
0 
0 


fa + Xaift ® ® 

0 fa + Aaafi $» + Asaft 
0 0 0 


0 0 

0 0 

fi + A.ift ft + *uft 


or 


8 . 

'(98) 

X = A x g, * = TA x g. (") 

Since T is constant we can do Step 3 of Table 3 directly on the natural curvatures. 
Now Ar(fiift.ft) i* » Kww function of the triangular coordinates. Consequently, the 

£!«..«. «- b. -imply «bUi»«J <2 ““ '“i /* ", 

- fy m 1/3. Let the corresponding matrix be A x . Then X = A x g, and the natural 
deviatoric curvatures are given by 

Xj = (a,-T x )», < 100 > 

which transformed to deviatoric Cartesian curvatures =* * - * gives finally 

= T(Av - A x )g=* A^g- ( 101 ) 


11.6 The Elem ent Stiffness Matrix 

The bask stiffness matrix K„ is the same derived in Ref. 8 using the conventional 
FF and need not be rederived here. The higher order stiffness matrix is given by Eqs. 84, 
which for a plate bending element specialise to 

K„ = aQ r KodQ * «Q T [f A AjDA^dA] Q, (102) 

where D is the Cartesian moment-curvature constitutive matrix resulting from the inte- 
gr&tion of E through the plate thickness: 

Du -^13 Dial f 
Dia Dyx Dia I < *** 

D\3 D 23 D $ 3 J \ v 

Since Ad varies linearly, if D is constant we could numerically integrate in 
Eq. 102 exactly with a three point Gauss rule, for example the three midpoint formula. 
But as the formation of the element stiffness is dominated by these calculations it is of 
interest to derive in closed form. Such a derivation is provided in Ref. 17. 

11.7 Preliminary Evaluation 

A. of this writing only a sketchy evaluation of the first ANDES element is available. 
We have found that for triangle, with good aspect ratio its behavior is similar to that of 
the scaled FF element of Ref. 8, which is known to be an excellent performer. But the 
ANDES element shows less distortion sensitivity for high aspect ratio elements, as can be 
expected from its construction. Additional evaluation details will be reported in Ref. 17. 

These preliminary results are encouraging in that we now have two good stand-alone 
components (FF and ANDES) of K*. Thus it is plausible that a weighted mix of these 
formulations as per Eq. 72 can be used to squeese the ultimate in performance for this very 
simple element. 
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-12. CONCLUSIONS 


1 . 

2 . 

3. 

4. 

5 . 


6. 

7. 


The results presented in this paper may be summarised a. follows. 

The classical variational principles of linear elasticity may be embedded m a 

T^^Usridty principles with independently varied displacements are members of 
a three-parameter family. Those principles without independent duplacementa are 

b * 

W n«ntlv investigated on this family nsmg hybnd forcing potentials. 

Kiwwti, .on.tt.mU ..Lid. <1~ « «f *• J ”** 

SETi— . _ U 

tbi.nfm.nttd virUtional tttttaf. A •»•“? «* ANS, •=•“ ANDM, «t* ip* 
uraSy the decomposition of the stiffness equations into basic and higher order parts. 
In addition, combined FF/ ANDES forms emerge from the general parametrised 

Thesatbfaction of the individual element test yield, various orthogonality conditions 
that the kinematic constraints should satisfy a prion. 

The first ANDES element based on this formulation displays encouraging s an - 
i regard, distortion sensitivity. The weighted combination of 
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